Đo lường tự tương quan trong R

Ma trận láng giềng không gian

Theo INSEE trang 48

sử dụng gói spdep

Cách mới nhất là theo Paulamoraga

Thư viện

#library(spData)
library(sf)
library(spdep)
library(ggplot2)
bando <- st_read("F:/five/R/thuchanhR/shape/bacninh.shp")

Reading layer `bacninh' from data source `F:\five\R\thuchanhR\shape\bacninh.shp' using driver `ESRI Shapefile'
Simple feature collection with 8 features and 8 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 105.9025 ymin: 20.96912 xmax: 106.3119 ymax: 21.2643
Geodetic CRS:  WGS 84

Xem bản đồ

plot(bando[2])

Định nghĩa láng giềng dựa trên tiếp giáp

library(spdep)
nb <- spdep::poly2nb(bando, queen = TRUE)
head(nb)

[[1]]
[1] 4 6 8

[[2]]
[1] 3 4 5

[[3]]
[1] 2 5

[[4]]
[1] 1 2 5 6

[[5]]
[1] 2 3 4 6

[[6]]
[1] 1 4 5 7 8

Vẽ bản đồ

plot(st_geometry(bando), border = "lightgray")
plot.nb(nb, st_geometry(bando), add = TRUE)

Warning in st_point_on_surface.sfc(coords): st_point_on_surface may not give
correct results for longitude/latitude data

Tiếp giáp với đơn vị diện tích số 5, id =5

id <- 5 # area id
bando$neighbors <- "other"
bando$neighbors[id] <- "area"
bando$neighbors[nb[[id]]] <- "neighbors"
ggplot(bando) + geom_sf(aes(fill = neighbors)) + theme_bw() +
  scale_fill_manual(values = c("gray30", "gray", "white"))

Liên kết dữ liệu điều tra vào shape file ?

Tải dữ file điểm nhà

Xem bản đồ

Láng giềng dựa trên k láng giềng gần nhất

Thư viện thống kê

# Neighbors based on 3 nearest neighbors, càng lớn càng có nhiều láng giềng
coo <- st_centroid(bando)

Warning: st_centroid assumes attributes are constant over geometries

nb <- knn2nb(knearneigh(coo, k = 2)) # k number nearest neighbors
plot(st_geometry(bando), border = "lightgray")
plot.nb(nb, st_geometry(bando), add = TRUE)

Warning in st_point_on_surface.sfc(coords): st_point_on_surface may not give
correct results for longitude/latitude data

Láng giềng dựa trên khoảng cách

# Neighbors based on distance: trong đường kính d1 đến d2 Từ huyện này sang kia tối thiểu 7km, xa nhất 12
nb <- dnearneigh(x = st_centroid(bando), d1 = 6, d2 = 12)

Warning: st_centroid assumes attributes are constant over geometries

plot(st_geometry(bando), border = "lightgray")
plot.nb(nb, st_geometry(bando), add = TRUE)

Warning in st_point_on_surface.sfc(coords): st_point_on_surface may not give
correct results for longitude/latitude data

Tìm sao cho mỗi khu vực có ít nhất 1 hàng xóm qua 2 bước

# Lấy trọng tâm, chuyển  đổi danh sách thành lớp với các id láng giềng
coo <- st_centroid(bando)

Warning: st_centroid assumes attributes are constant over geometries

# k is the number nearest neighbors
nb1 <- knn2nb(knearneigh(coo, k = 1))

Warning in knn2nb(knearneigh(coo, k = 1)): neighbour object has 2 sub-graphs

# Sau đó truyền danh sách nb và tính toán tóm tắt
dist1 <- nbdists(nb1, coo)
summary(unlist(dist1))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  6.130   7.165   7.610   7.732   8.250   9.733 

# Kết quả d1=6.130 và d2 = 9.733

Láng giềng của bậc k dựa trên tiếp giáp

library(spdep)
nb <- poly2nb(bando, queen = TRUE)
nblags <- spdep::nblag(neighbours = nb, maxlag = 2) # tạo danh sách chứa láng giềng bậc 1 và bậc 2

# Neighbors of first order
plot(st_geometry(bando), border = "lightgray")
plot.nb(nblags[[1]], st_geometry(bando), add = TRUE)

Warning in st_point_on_surface.sfc(coords): st_point_on_surface may not give
correct results for longitude/latitude data

Láng giềng bậc 2

# Neighbors of second order
plot(st_geometry(bando), border = "lightgray")
plot.nb(nblags[[2]], st_geometry(bando), add = TRUE)

Warning in st_point_on_surface.sfc(coords): st_point_on_surface may not give
correct results for longitude/latitude data

Hàm spdep::nblag_cumul() tích lũy một danh sách láng giềng bị trễ và nblag() trả về danh sách láng giềng duy nhất.

# Neighbors of order 1 until order 2
nb <- spdep::poly2nb(bando, queen = TRUE)
nblagsc <- spdep::nblag_cumul(nblags)
plot(st_geometry(bando), border = "lightgray")
plot.nb(nblagsc, st_geometry(bando), add = TRUE)

Warning in st_point_on_surface.sfc(coords): st_point_on_surface may not give
correct results for longitude/latitude data

Ma trận láng giềng

Ma trận trọng số không gian dựa trên danh sách láng giềng dạng nhị phân

Có nhiều dạng style = B là mã hóa nhị phân cơ bản và W là được chuẩn hóa theo hàng (1/số láng giềng)

nb <- poly2nb(bando, queen = TRUE)
nbw <- spdep::nb2listw(nb, style = "W")
nbw$weights[1:3]

[[1]]
[1] 0.3333333 0.3333333 0.3333333

[[2]]
[1] 0.3333333 0.3333333 0.3333333

[[3]]
[1] 0.5 0.5

Hình dung ma trận trọng số không gian bằng cách tạo ma trận có trọng số bằng listw2mat() và sử dụng lattice::levelplot() để vẽ

library(lattice)
m1 <- listw2mat(nbw)
lattice::levelplot(t(m1),
scales = list(y = list(at = c(10, 20, 30, 40),
                       labels = c(10, 20, 30, 40))))

Ma trận trọng số không gian dựa trên giá trị nghịch đảo khoảng cách

coo <- st_centroid(bando)

Warning: st_centroid assumes attributes are constant over geometries

nb <- poly2nb(bando, queen = TRUE)
dists <- nbdists(nb, coo)
ids <- lapply(dists, function(x){1/x})

nbw <- nb2listw(nb, glist = ids, style = "B")
nbw$weights[1:3]

[[1]]
[1] 0.08803397 0.13315812 0.08922956

[[2]]
[1] 0.16312086 0.12604479 0.06963587

[[3]]
[1] 0.16312086 0.06395368

m2 <- listw2mat(nbw)
lattice::levelplot(t(m2),
scales = list(y = list(at = c(10, 20, 30, 40),
                       labels = c(10, 20, 30, 40))))

Tự tương quan không gian

Xem mẫu

library(spData)
library(sf)
#library(mapview)
library(tmap)

Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
remotes::install_github('r-tmap/tmap')

#map <- st_read(system.file("shapes/boston_tracts.shp",
#                           package = "spData"), quiet = TRUE)
map <- st_read("F:/five/R/thuchanhR/shape/bacninh3.shp")

Reading layer `bacninh3' from data source `F:\five\R\thuchanhR\shape\bacninh3.shp' using driver `ESRI Shapefile'
Simple feature collection with 125 features and 13 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 105.9034 ymin: 20.96926 xmax: 106.3112 ymax: 21.26444
Geodetic CRS:  WGS 84

map$vble <- rnorm(125,20,5)
#map$vble <- map$MEDV # Giá trị trung bình bán nhà trong khu vực
tm_shape(map)+ tm_polygons("vble")

Chỉ số toàn cục Moran’s I

Hàm kiểm định moran.test()

H_0: I \leq E[I] (tự tương quan âm hoặc không tự tương quan) H_1: I > E[I] (tự tương quan dương)

# Neighbors
library(spdep)
nb <- poly2nb(map, queen = TRUE) # queen shares point or border
nbw <- nb2listw(nb, style = "W")

# Global Moran's I
gmoran <- moran.test(map$vble, nbw,
                     alternative = "greater")
gmoran

    Moran I test under randomisation

data:  map$vble  
weights: nbw    

Moran I statistic standard deviate = -0.39394, p-value = 0.6532
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
     -0.029881112      -0.008064516       0.003067050 

Mô phỏng

Đồ thị scatter về Moran I

moran.plot(map$vble, nbw)

Chỉ số Moran I cục bộ

lmoran <- localmoran(map$vble, nbw, alternative = "greater")
head(lmoran)

            Ii         E.Ii      Var.Ii       Z.Ii Pr(z > E(Ii))
1  0.765817899 -0.023805987 0.952557826  0.8090486     0.2092436
2  0.009519732 -0.000171506 0.002912713  0.1795686     0.4287456
3  0.149080936 -0.004442026 0.065165918  0.6013995     0.2737870
4  0.192576711 -0.008179690 0.196223798  0.4532038     0.3252010
5  0.351715687 -0.003373169 0.102493614  1.1091458     0.1336836
6 -0.503870117 -0.001514096 0.061967536 -2.0180387     0.9782064

Các thống kê về Moran I cục bộ

library(tmap)
tmap_mode("view")

tmap mode set to interactive viewing

map$lmI <- lmoran[, "Ii"] # local Moran's I
map$lmZ <- lmoran[, "Z.Ii"] # z-scores
# p-values corresponding to alternative greater
map$lmp <- lmoran[, "Pr(z > E(Ii))"]

p1 <- tm_shape(map) +
  tm_polygons(col = "vble", title = "vble", style = "quantile") +
  tm_layout(legend.outside = TRUE)

p2 <- tm_shape(map) +
  tm_polygons(col = "lmI", title = "Local Moran's I",
              style = "quantile") +
  tm_layout(legend.outside = TRUE)

p3 <- tm_shape(map) +
  tm_polygons(col = "lmZ", title = "Z-score",
              breaks = c(-Inf, 1.65, Inf)) +
  tm_layout(legend.outside = TRUE)

p4 <- tm_shape(map) +
  tm_polygons(col = "lmp", title = "p-value",
              breaks = c(-Inf, 0.05, Inf)) +
  tm_layout(legend.outside = TRUE)

tmap_arrange(p1, p2, p3, p4)

Variable(s) "lmI" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.

Variable(s) "lmZ" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.

Variable(s) "lmI" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.

Variable(s) "lmZ" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.

Tương quan cục bộ theo vùng trên bản đồ

tm_shape(map) + tm_polygons(col = "lmZ",
title = "Local Moran's I", style = "fixed",
breaks = c(-Inf, -1.96, 1.96, Inf),
labels = c("Negative SAC", "No SAC", "Positive SAC"),
palette =  c("blue", "white", "red")) +
tm_layout(legend.outside = TRUE)

Các cụm

Cao - Cao:

Cao - Thấp:

Thấp - Cao:

Thấp - Thấp:

lmoran <- localmoran(map$vble, nbw, alternative = "two.sided")
head(lmoran)

            Ii         E.Ii      Var.Ii       Z.Ii Pr(z != E(Ii))
1  0.765817899 -0.023805987 0.952557826  0.8090486     0.41848720
2  0.009519732 -0.000171506 0.002912713  0.1795686     0.85749125
3  0.149080936 -0.004442026 0.065165918  0.6013995     0.54757392
4  0.192576711 -0.008179690 0.196223798  0.4532038     0.65040198
5  0.351715687 -0.003373169 0.102493614  1.1091458     0.26736729
6 -0.503870117 -0.001514096 0.061967536 -2.0180387     0.04358723

Lược đồ

map$lmp <- lmoran[, 5] # p-values are in column 5

mp <- moran.plot(as.vector(scale(map$vble)), nbw)

head(mp)

           x          wx is_inf labels      dfb.1_       dfb.x       dffit
1 -1.7112355 -0.44394320  FALSE      1 -0.11598463  0.19927572 -0.23057157
2 -0.1452467 -0.06501746  FALSE      2 -0.02293651  0.00334486 -0.02317912
3  0.7391919  0.20006751  FALSE      3  0.03859288  0.02864234  0.04806032
4  1.0030789  0.19044972  FALSE      4  0.03837016  0.03864318  0.05445699
5 -0.6441481 -0.54164869  FALSE      5 -0.12779197  0.08264821 -0.15218908
6  0.4315621 -1.15820907   TRUE      6 -0.25791745 -0.11175532 -0.28108835
      cov.r       cook.d         hat
1 1.0221767 0.0264465120 0.031615539
2 1.0237387 0.0002706931 0.008170134
3 1.0261344 0.0011626113 0.012406489
4 1.0300486 0.0014927206 0.016114252
5 0.9949362 0.0114856797 0.011346183
6 0.9005201 0.0373103399 0.009501983

Tạo biến quadrant thể hiện 4 góc lược đồ quan hệ tương quan cục bộ

map$quadrant <- NA
# high-high
map[(mp$x >= 0 & mp$wx >= 0) & (map$lmp <= 0.05), "quadrant"]<- 1
# low-low
map[(mp$x <= 0 & mp$wx <= 0) & (map$lmp <= 0.05), "quadrant"]<- 2
# high-low
map[(mp$x >= 0 & mp$wx <= 0) & (map$lmp <= 0.05), "quadrant"]<- 3
# low-high
map[(mp$x <= 0 & mp$wx >= 0) & (map$lmp <= 0.05), "quadrant"]<- 4
# non-significant
map[(map$lmp > 0.05), "quadrant"] <- 5

Thể hiện trên bản đồ

tm_shape(map) + tm_fill(col = "quadrant", title = "",
breaks = c(1, 2, 3, 4, 5, 6),
palette =  c("red", "blue", "lightpink", "skyblue2", "white"),
labels = c("High-High", "Low-Low", "High-Low",
           "Low-High", "Non-significant")) +
tm_legend(text.size = 1)  + tm_borders(alpha = 0.5) +
tm_layout(frame = FALSE,  title = "Clusters")  +
tm_layout(legend.outside = TRUE)